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Abstract 

The Ruderman-Kittel-Kasuya-Yosida (RKKY) interaction between two 
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\Q ' magnetic s-d spin impurities across a tunneling junction is studied when the 
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system is driven out of equilibrium through biasing the junction. The nonequi- 
librium situation is handled with the Keldysh time-loop perturbation formal- 



ism in conjunction with appropriate coupling methods for tunneling systems 

o" 

J-^ ' due to Caroli and Feuchtwang. We find that the presence of a nonequilib- 

rium bias across the junction leads to an interference of several fundamental 
oscillations, such that in this tunneling geometry, it is possible to tune the 
interaction between ferromagnetic and antiferromagnetic coupling at a fixed 
impurity configuration, simply by changing the bias across the junction. Fur- 
thermore, it is shown that the range of the RKKY interaction is altered out of 
equilibrium, such that in particular the interaction energy between two slabs 
of spins scales extensively with the thickness of the slabs in the presence of 
an applied bias. 

PACS numbers: 75.30Et, 73.40Gk, 85.30Mn, 75.70-i 
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I. INTRODUCTION 



The Ruderman-Kittel-Kasuya-Yosida (RKKY) interaction has been a very intensely 
studied phenomenon in solid state physics since it was first proposed as an interaction 
between nuclear spins0i, and later between localized electronic spins in metalsi. More re- 
cently studies of the interaction have focussed on its effects in various structured systems 
and in particular on the role it plays in the giant magnetoresistance effects observed in some 
layered structures of magnetic and non-magnetic materials^i. 

Furthermore it has been suggested that the RKKY interaction is responsible for spin po- 
larization effects observed in tunneling systems!, such as layered structures with a potential 
barrier formed by an insulating oxide layer or a vacuum gap. For arrangements involv- 
ing movable (vacuum) tunneling junctions, such as those occurring in scanning tunneling 
microscopy (STM) it has been suggested that the exchange interaction between two mag- 
netic materials on either side of the tunneling junction can be used in a modified version of 
atomic force microscopy, which may be called exchange force microscopy, in order to resolve 
an atomic image of a sample structure&0. 

Particularly in tunneling systems, the measurement of electronic properties, such as for 
example the electronic density of states, intrinsically requires the structure to be biased 
out of equilibrium. The presence of a nonequilibrium bias may also occur in structures 
which exhibit giant magnetoresistance phenomena, when parts of the structure contain a 
potential barrier. To date, however, descriptions of the RKKY interaction are confined 
to systems in equilibrium, including approximate theoretical treatments of the interaction 
across a potential welfl and a tunneling barrierS. 

It is therefore the purpose of the present paper to establish a theoretical description of 
the RKKY-interaction across a tunneling barrier out of equilibrium from first principles. For 
this purpose we employ the Keldysh nonequilibrium perturbation formalismO along with a 
coupling procedure for structured systems developed by Caroli, Combescot et a/.0Hll and 
FeuchtwangEUll, leading to a proper nonequilibrium field theoretic description. Since the 
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treatments by Caroli and Feuchtwang are proper many-body formalisms their application to 
the present problem simultaneously provides the basis for the inclusion of further many-body 
effects such as carrier-carrier or carrier-phonon interactions to the problem. 

With this technique we manage to derive a general nonequilibrium solution for the RKKY 
interaction for various dimensionalities and arrangements of spins within a structured sys- 
tem. Particularly we show that the interaction of spins across a tunneling junction can be 
tuned between ferromagnetic (FM) and antiferromagnetic (AFM) coupling by changing the 
bias across the junction alone. For the application to exchange force microscopy, mentioned 
before, it is shown that varying the bias across the vacuum junction in an STM can lead to 
a force which switches between attractive and repulsive behavior and that this force should 
be experimentally measurable with a state of the art apparatus. Furthermore, we show that 
the presence of a nonequilibrium bias across the junction significantly alters the range of 
the RKKY interaction, such that in particular the interaction energy between two slabs of 
spins scales extensively with the thickness of the slabs. 

The remainder of the paper is structured as follows: In Sec. II we establish a model of 
the tunneling system containing the spin impurities. Sec. Ill contains a derivation of the 
RKKY interaction out of equilibrium in various dimensions in terms of general Keldysh 
nonequilibrium Green's functions. In Sec. Ill A we consider the interaction of two magnetic 
s-d impurities across a one dimensional tunneling junction out of equilibrium and in Sec. 
Ill B we extend these results to the more realistic situation of two magnetic layers or 
slabs of spins on either side of a planar three dimensional tunneling junction, including 
the possibility of a non-magnetic spacer material between the magnetic materials and the 
barrier. Sec. IV describes how the Keldysh Green's functions used in Sec. Ill are obtained 
in a system containing a potential barrier. In Sec. V we implement our results numerically 
and calculate equilibrium and nonequilibrium versions of the RKKY interaction for various 
tunneling geometries. Sec. VI discusses the experimental observability of the behavior of the 
interaction predicted by the numerical study. In Sec. VII, in conclusion, the implications 
and possible further applications of this work are summarized. 
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II. THE MODEL 



The one dimensional tunneling system we consider consists of a tunneling barrier which 
is connected to two leads on the left and right at the points L and R, respectively, as shown 
in Fig. 1. In equilibrium the barrier is assumed to be flat on top with an abrupt potential 
change of V at the interfaces. The corresponding single particle potential can be written as 

V(x) = V Q(x-L)e(R-x). (1) 

Upon biasing the junction, the potential of the barrier acquires a slope and the chemical 
potential fi R in the right lead undergoes a shift eV with respect to the chemical potential 
fj, L in the left lead. Simultaneously, the conduction band bottom V R in the right lead shifts 
by eV with respect to conduction band bottom V L in the left lead, such that the electronic 
potential is changed to 

V(x) = \v - eV^—^\ B(x - L)Q(R - x) - eVQ(x - R). (2) 
R — Li 

Here the bottom of the conduction band on the left side is taken as the origin of energy, e 
is the modulus of the elementary charge of an electron, and V the voltage drop across the 
barrier. 

Two magnetic s-d impurities are situated within the electrodes, on either side of the 
barrier, at an equal distance d from the interfaces. For simplicity, both electrodes are 
considered to consist of the same material and the effective electron masses are assumed to 
be equal in all three parts of the junction. 

The Hamiltonian of the system, with a single electron band in each part of the junction, 
can be written as 



H = £ / dx¥ a (x) 



— + V(x) 



E CT «/3 • S p / dx¥ a (x)5(x - x p )V p {x) =H + H s _ d , (3) 
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p=l,2;a,/3 



where H s _d is the last, spin dependent, term in the Hamiltonian and the coupling constant 
J is assumed to have the units Jm d , where d is the dimensionality of the system considered. 



Here a space representation has been chosen, which proves to be advantageous due to the lack 
of translational invariance caused by the presence of the barrier potential V(x). p 2 /2m is the 
kinetic energy of electrons with uniform effective mass m, and we use h = 1 throughout. (T a p 
is the vector of Pauli matrices, used to represent the spin of the conduction electrons which 
couple to the two local moments S p . The indices a and (3 label the two spin components: 

a,/3e{U}. 

III. EXPRESSION FOR THE RKKY INTERACTION 

A. The interaction in one dimension 

The RKKY interaction energy between the two impurities is calculated from the lowest 
order exchange contributions to the perturbation of the energy of the localized spins in the 
presence of the conduction electrons. Either in or out of equilibrium this energy is given 
by the expectation value (if s _d)(i ex ), where the subscript (lex) indicates that only the first 
order exchange contributions to the average are considered. Effectively, therefore, one has 
to calculate the first order perturbation to the conduction electron spin density due to the 
first spin at the location of the second spin and vice vers$$^. The total interaction energy 
is a sum of these two contributions 

^RK = -- S P ■ A P(lex)(z P ), (4) 

Z p=l,2 

where 

a,/3 a,0 

is the first order exchange contribution to the spin density at site p. 

In the present case we now have to calculate (pj) by means of ([5]) out of equilibrium. In 
order to do so we make use of the Keldysh formalism. In the Keldysh notation the spin 
dependent particle density correlation function (n a p(x)) = (\I / ^(x)\I' ( g(x)) can be written in 
terms of the Keldysh Green's function G^Jx, t; x', 0) which is defined as 
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G^(x,t;x',0) =i(*J(x , ,0)* a (x,t)) > (6) 

where the operators ^ a (x,t) and \1/ Q ,(x, t) are the field operators of the system considered, 
which create and destroy a particle with spin a at point x and time t, respectively. From 
this we find 

/no rlii) 
—e^G< ap {x,x'-u). (7) 

X 1 — HE 

Using the definition of G K , it can be shown in equilibrium that this expression is equivalent 
to the usual relation between the particle density and the retarded Green's function G r , 

1 f°° r ] 

(n a p(x)) — lim / n^(ci;)Im G r a Jx, x ; u>) du, (8) 



IT x'—*x 
-1 : 



where np(uj) = {exp[/3(tu — /x)] + 1} is the usual Fermi distribution function in equilibrium, 
fi is the equilibrium chemical potential and G r is the retarded Green's function of the system, 
which is defined as 

G r a p(x, x'- 1) = -iQ(t) 0), V a (x, t)}), (9) 

where { , } denotes the anti-commutator. The advanced Green's function G a 

G a a/3 (x, x'- 1) = ie(-t) ({^(x', 0), V a (x, t)}), (10) 

will be used later. 

In order to obtain a proper nonequilibrium result for the RKKY interaction we have to 
perform first order nonequilibrium perturbation theory on expression ([J]). The appropriate 
formalism, introduced by Keldysh, reformulates the regular diagrammatic perturbation the- 
ory in terms of a 2 x 2 matrix formalism, to properly handle the time development, since out 
of equilibrium the usual S'-matrix expansion based on the Gell-Mann-Low theorem breaks 
downffl. The Green's function matrix of the unperturbed system without spin impurities, 
but including the full barrier potential within this formalism is written as 



G(o)a/?(£, x'; t) 



y Gfo) a p{ x i X ''i _ ^(0)a/3( X ' X> 'l ) 
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(11) 



Cr* and Cr* denote the usual time-ordered and anti-time-ordered Green's functions with the 
general definitions 

G t a p{x, x'; t) = -i(TV a (x,t)^(x',0)), (12) 
G^{x, x'; t) = -i{f* a {x, t^lix', 0)>, (13) 

where T and T are the time-ordering and anti-time-ordering operators, respectively. G > is 
the complementary correlation function to G < : 

GZ p (x,x';t) = -i(* a (x,t)¥p(x',0)). (14) 

The corresponding matrix perturbation expansion to first order in the coupling J yields 

Gr(l) Q /3(Xp, x p ; t — t') = G( ) ai g(Xp, x p ; t — t')5 a ^ 

— — S q ■ cr p a I dtiQ(p) aa {x pi x q \t — ti)Q(Q)j3p{x qi Xp\ti — t'). (15) 

1 g=l,2 J 

From ( P^D the first order exchange terms of the function G < can be resolved as 

^\\cY^afi^ X Vi X Pi t ~~ = _ ~2^* q ' 17 P a J ^ l ^(0) ^ _ ^l)^(0) (^P) ^1 — 

+ G( 0) (x p ,a; o ;t - tx)Gf Q) (x p ,x q ',t\ - t')\ (pe{l,2},g^p), (16) 

where we have now dropped the spin indices in the unperturbed Green's functions since 
they are spin independent. To obtain (|16|), we have used appropriate relations between the 
six Green's functions G < , CP, Cr*, Cr*, G r and CP in order to replace the dependence of the 
result on G l and G l by one on G r and G a (cf. Ref. H). 



One relation arising from this transformation that holds generally for operators in the 
Keldysh formalism and which is particularly useful to us is 

(AB)< = A<B a + A r B<, (17) 

where AB is to be understood as a matrix product of two general operators A and B, imply- 
ing integrations over space and time where applicable. By means of Fourier transformation 
of (|16|) into the frequency domain, (f|) can therefore be expressed as 



E 



RK 



J 2 Si • s 2 



Im 



G(p)(xi, X2) 



x {G r {0) (x 2 ,xx) + G? 0) (x 2 ,x 1 )} 



(0)1 



duo 
2n' 



(18) 



The above equation is the general expression for the RKKY interaction out of equilibrium 
in a purely ID system without translational invariance. It is one of the central results of the 
present work and the principal goal of the following treatment is to evaluate it explicitly for 
various situations. 

In an equilibrium situation it is easily established that QIBp goes over to the well known 



resu 



E KK = — J 2 Si ■ S 2 / n F [ 

71 J-00 



x Im 



(0)> 



dw. 



(19) 



Before we establish how the Keldysh Green's functions Gf^ and G r J^ are calculated in a 



nonequilibrium situation we will generalize the results just obtained to higher dimensions. 



B. Generalization to two and three dimensions 

For the extension to higher dimensions we assume that the system is translationally 
invariant in the further one or two dimensions, corresponding to a perfectly planar barrier. 
Once we have found the solution for the purely ID Green's function G^, we can simplify the 
solution of the planar extension of the equivalent problem by exploiting the translational 
invariance of the system in the additional directions parallel to the barrier by means of 
corresponding Fourier transforms 

G^V^kii;^) = J d rf - 1 re lk n r G 2 / ) 3D (x,x», (20) 

where x is a d-dimensional coordinate, r = xy — xji is the (d— 1) -dimensional relative position 
vector parallel to the barrier and kn denotes the (d — l)-dimensional electron wavevector 
parallel to the barrier, where d £ {2, 3}. In the following we shall drop the superscripts again 
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that were just introduced to indicate the dimensionality of the system which the correspond- 
ing Green's functions describe, since it will be always recognizable from the arguments of 
these functions whether they pertain to a purely ID system or to a planar version in higher 
dimensions. 

The extension to 2/3D leaves the spatial dependence of the Hamiltonian unaffected in 
the direction perpendicular to the barrier and, as will be shown in Sec. IV, only amounts to 



a change in the energy arguments of the corresponding Green's functions (see also Ref. |18|) . 



The expression for the RKKY interaction between two arbitrarily placed impurities with 



respect to the barrier in (i-dimensions can then be written as an extension from ([18 



h RK - J bx ■ b 2 J {27i)i 2d-2) e 11 



(2tt)( m - 2 ) 

/oo r 
Im Gf 0) (x u x 2 ;q\\ - k^u) 
-oo 

du 

(dG{2,3}), 



x {G[ 0) (x 2 ,xi;k||;u;) + G a {0) (x 2 ,x 1 ;k ll ;uj)}] — (21) 



where r 12 denotes the impurity displacement in the direction parallel to the barrier. 

The most immediate application of the present theory is in a 3D tunneling system where 
the interaction between either two monolayers or two slabs of spins across a barrier is of 
interest for the description of giant magnetoresistance phenomen It will be assumed 
that all spins within one monolayer or slab have the same orientation due to a predominant 
ferromagnetic coupling over short distances. In equilibrium and without a barrier it was 
pointed out by Yafeti! how the interaction between two monolayers can be obtained from 
the conventional 3D version. In our case the corresponding expression is obtained through 
integrating over all r i2 and subsequently over qy in (|2l|) : 

/ \ 2 r d^kw r°° r 

£rk= (jp 2 s D d ) ■ S 2 J -J- 2 J Jm [Gf 0) (x u x 2 ; -k„; W ) 



x {^(ojC^^ijkiijw) + G^ 0) (x 2 ,x 1 ;k||;a;)}l (22) 



(27T) 

2tt' 

where E^ K is now the monolayer interaction energy across the junction, v 2D is the surface 
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area of the junction and p 2 s _ d is the surface density of the s-d spins in each of the two planes 
parallel to the barrier. 

The situation where two finite slabs of spins interact across a finite spacer layer which 
contains the barrier is shown schematically in Fig. The interaction between two slabs of 
spins across a planar tunneling barrier can be obtained by summing the contributions coming 
from each pair of interacting monolayers involved. In the present case where a continuous 
system is considered this amounts to two spatial integrations of the monolayer interaction 
in (P^ ) in the direction perpendicular to the barrier. If, for simplicity, we consider two 
interacting slabs of the same thickness I which are symmetrically displaced from the barrier, 
these two integrations assume the form 

/ *n \ 2 f L ~ d fR+(d+l) 

Kk = [Ps-dJ , d Xl dx 2 

v ' JL-{d+l) JR+d 

where E^ K is the interaction energy of the two slabs and in E^ K (xi,X2), is taken from ( |22|) 
where x\ indicates a position within the left electrode and x 2 one within the right one, 
respectively, and p^ d is the 3D density of s-d spins within the slabs. 

Now that formal expressions for the RKKY interaction are in hand, the next step is the 
evaluation of the Keldysh Green's functions Gf Jx, x') and G r ^(x, x') of the nonequilibrium 
ID system and their extensions to higher dimensions. 



E% K (xx,x 2 ) 



2D 
Ps-d 



(23) 



IV. EVALUATION OF THE RELEVANT GREEN'S FUNCTIONS 

It has been established by Carol&il and FeuchtwangEHl, that a tunneling system out 
of equilibrium can be treated by partitioning the system into several uncoupled parts which 
are considered to be in equilibrium at t = — oo. These parts are subsequently coupled to 
each other through appropriate transfer terms, in conjunction with an adiabatic switching 
procedure, to finally yield a nonequilibrium steady state. Formally, this corresponds to a 
perturbation expansion of these transfer terms to all orders@lll, but it can also be shown 
to be equivalent to applying Green's theorem at the partitioning points§'0. The number of 
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partitions to be made in the system is in principle arbitrary and will largely depend on the 
geometry considered. In the present case with a possibly sloping tunneling barrier of finite 
width, a partitioning into three regions at the electrode-barrier interfaces L and R - shown 
by the vertical dashed lines in Fig. 1 - is most convenient. In the present treatment we will 
follow closely the approach of Ref. [T7| for continuous systems. The Hamiltonian H in (|3|) is 
consequently written as a sum of three independent parts 

H (x) = G(L - x)H L {x) (24) 
+ 0(x - L)Q(R - x)H B {x) + 0(x - R)H R (x). 

A. Green's functions in one dimension 

In one dimension the Green's function G( ) of the system including the barrier has to 
satisfy the inhomogeneous Schrodinger equation 

[uo - H (x)}G (0) (x,x';uj) =6(x-x'). (25) 

Similarly one can define Green's functions for the several uncoupled sub-parts of the system 
i] e {L,B,R} as 

[uj — H^x^g^x, x'; u) = 5(x — x'), (26) 

where x, x' lie within the appropriate region determined by the choice of rj. Additionally, 
these functions have to satisfy appropriate boundary conditions at the electrode-barrier 
interfaces R and L. In principle these boundary conditions can be an arbitrary mixture 
of the Dirichlet and Neumann type. However, for simplicity it is best to use one or the 
other exclusively, and we choose here the Dirichlet conditions that the Green's functions 
vanish if they are taken with one of their arguments on the respective interfaces R or L. 
Relations between these uncoupled Green's functions g v and the Green's function G( ) of the 
full system, satisfying (fZ5"D, can be established, as noted before, by means of Green's theorem 
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applied at the interfaces. If in addition the discontinuity conditions for the derivatives of 



the full Green's function Gt \ are used 



2m = d x >Gr )(x,x' 



x'=x~ 



(27) 



(where we have left out the energy argument for brevity), together with appropriately chosen 
versions of further continuity conditions 



G( )(x,x') 



x=x 
x=x' 



d x d x >G(p)(x,x') 



0. 



(28) 
(29) 



which this function has to satisfy, one can find G(p)(x,x'). In the simple cases where x, x' G 
{L, R}, the full Green's function can be written a 



( G {0) (L,L) G {0) (L,R)\ 



J 



( 



2m 



2m 
~D 



\G(p){R,L) G(p)(R,R) 

lL (L,L)+ lB (L,L) -y B (L,R) 
v -1b{R,L) 7r (R j R)+ 1b (R,R) J 

7 R (R,R)+7b(R,R) 1b{L,R) 

j b (R,L) Jl(L, L) + 7b (L, L) 



(30) 



where 



D = [ lR {R, R) + 1B (R, R)][ 1L (L, L) + lB (L, L)\ 
-j B (L,R)j B (R,L) 



(31) 



and 



j v (a,b) = -—d x d x >g v (x,x';u) \x=a . 

2m x'=b 



(32) 



The above Green's function is readily cast into a retarded or advanced version by analytically 
continuing uj — > lim,5^o+ uj ±i5. 

The calculation of G^, however, is significantly complicated through the matrix property 
(|17|) inherent to all its defining equations in terms of continuity conditions. The appropriate 
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choice for these continuity conditions are similar to the ones for G r ^ with the only important 
difference that all first derivatives are continuous also now at x — x' . For a general situation 
with two partitions we find that 



Gfo) (-k> L) G< j (L, R) ^ 



2m 



J 



D 



(33) 



x 



x 



X 



/ 



V 



\ 

/ 

\ 

) 

\ 

) 



yG< 0) (R,L) G< 0) (R,R) 
Y R (R, R) + Yb(R, R) Yb(l,R) 

Yb(R,L) Yl(L,L)+Y b (L,L) 

7 f(L,L)+7l(A^) ~lUL,R) 

-j§(R,L) 1 <(R,R)+ 1 <(R,R) 

Y R (R,R) + Y B (R,R) Yb(L,R) 

Yb(R,L) Yl(L,L) + Yb(L,L) 

where it should be noted that despite the presence of the terms with 7^ in the matrix 
in the middle, the full Green's function G^s does not depend on these terms. In a true 
tunneling situation, where one only considers energies lower than the height of the barrier, 
this is immediately seen to be true since for those energies no states exist in the barrier 
region. Another case, however, is the one where there are states inside the barrier region 
in the energy range considered, such as impurity levels or quasi-bound states in a double- 
barrier tunneling structure; or when a scattering problem across a non-periodic potential is 
considered rather than a tunneling situation. For these cases one can showQ that terms i 
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G^ which contain 7^, occur in conjunction with terms coming from |-D| and other terms 
from the matrix product in (|33|) such that they fulfill an identity A5(A)A = 0, where A is 
the denominator of 7^, which is independent of the spatial arguments of 7^. Put another 
way, the contributions of the poles in 7^ which constitute 7^ are suppressed in G^, since 
wherever they occur they are given zero weight. 

One can understand this cancellation from a physical picture of how the nonequilibrium 
system is established. While the left and right leads are modeled as semi-infinite grand 
canonical ensembles with (possibly different) chemical potentials on either side at t = — 00, 
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the barrier region is finite and not coupled to any exterior reservoir. Once a steady state 
is established, any quantity of the fully coupled system will depend only on the initial 
occupations of the semi-infinite leads. In particular, there can be no dependence of on 
7_b, since the infinite, fully coupled system cannot remember the initial occupation of the 
finite barrier region. In the same way it turns out that in terms containing r y r ^ a only the 
real part jb contributes to Gf^i as contributions coming from the imaginary parts cancel 
for the reasons outlined above. Therefore we can set 7^ = and leave away the superscripts 
in 7^ a for further calculations. 

For the evaluation of (0), however, we still need to know the appropriate expressions 
for 7<, 77 G {R, L}. Since the left and right decoupled regions are separately in equilibrium, 
they have well defined electron occupations n F (u), and therefore we find that g< can be 
expressed as 

g<(x, x'\ u) = n' F (uj) [g°(x, x'; w) - g*(x, x'; w)] , (34) 

which immediately transfers to the 7^ as 

7<(x,x';w) = n F {u) [rf(x,x';u) - rf(x, x'\ u)] . (35) 

If the impurities lie deeper within the electrodes on the left and right hand sides of 
the barrier, the relevant Green's functions G r J^(xx,x 2 ) and G^(xx,x 2 ) can be expressed in 
terms of the full Green's functions between the interfaces, (^) and (|33|), and the Green's 
functions of the uncoupled leads, 

G r ( ^(xi,x 2 ) = -(2m)~ 2 d xl g r £ a {x u x') ^ =l _ (36) 



r (o) 
r (o) 



x G t ,q?(L, R) d x ,g R /a (x',x 2 ) 



Gf 0) (x 1 ,x 2 ) = -(2m)- 2 (37) 



x {d x >gt{xi,x')\ x > =L - G a (0) {L,R) d x ,g a R (x',x 2 )\ x , =R+ 



r (o) 

r (0) 

+ d x/ g r L (x u x')\ x , =L . Gf 0) (L,R) d x ,g R (x\x 2 )\ x , =R+ 
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+ d x/ g r L (x u x')\ x , =L - G r (0) (L,R) d x >g^{x', x 2 )\ x , =R+ } 
(xi < L, X2 > R). 

The corresponding expressions for the Green's functions with reversed arguments are ob- 
tained in an analogous way. 

The RKKY interaction (pD can now be expressed entirely in terms of the unperturbed 



Green's functions of the separate subsystems in equilibrium using (^) and (|33|)-([37l). For 
the simple case where the impurities are situated immediately on the left and right barrier- 



electrode interfaces (18) can for example be expressed in terms of the quantities / y v as 

E RK = -2(2m) 2 J 2 Sx ■ S 2 J ^|D| _2 Re[D _1 ]7B(L, R)jb(R, L) 
x Im{[7*CR, R) + MR, R)} L) + j£(R, R) [ML, L) + Jl(L, L))} . (38) 

Within the single effective mass approximation for a barrier system as shown in Fig. [I] 
the functions of the decoupled leads are found to be 

g$n i M = +(-)^L (39) 

1l(r) 

sin[g^ ) (x-L( J R))]e T(±) ^) (x '~ L(jR)) x> {<)x' 



x 



sm[q r ^ a R) (x' - L(i?))]e T(±) ^(«) (x ~ L ™ x'> {<)x 



where ?2(fl) = ^2m(w — V L ^ ± i5), the upper/lower signs are associated with the super- 
scripts r/a, respectively, and V L<yFC> is the bottom of the conduction band in the corresponding 
side of the junction as before. Likewise the corresponding Green's functions for the sloping 
barrier region are obtained as 

, 2m7r/t~ 1 
SB(X - X)= M)HL) - f(L)h(R) (40) 



X 



[h{L)f(x) - h{x)f(L)} [f(R)h{x') - f{x')h(R)} x<x' 
[h(R)f(x) - h(x)f(R)} [f(L)h(x>) - f{x')h{L)\ x>x> 



where 
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f(x) = Ai(Kx + (/ri 2 ), 
h(x) = Bi(«aj + C/v 2 ), 

are two independent solutions of the inhomogeneous Schrodinger equation 



(41) 



dl - k 3 x 



f{x) = C, 



(42) 



3 /2meV and £ = 2m {y Q _ e yi 2 _ u y 



with the parameters k — v 

In equilibrium, where there is no slope to the barrier, one can show that qb{x,x') sim- 
plifies to 

sinh[fc(x — L)] smh[k(x' — R)] , 

(43) 



9b(.x,x 



2m ( ks\v&[k(R- L)\ 

sinh[fc(x' — L)\ sinh[/c(x — R)\ 



ksmh[k{R-L)] 



x' > x 



where k = y2m[Vo — lu]. Note that (|40|) and (|43|) hold for all oo, i.e. including the case when 
Vq — uj < 0, for which the sinh-functions in (^) go over to corresponding sin-functions. 

When the results of the present and the previous section are combined in equilibrium 
and for the limit of the barrier height or the barrier width going to zero we can obtain 
analytic results for all versions of the interaction considered so far. The purely ID result 
from ( |l~8t ) for this case reduces to the well known expression for the RKKY interaction in 
one dimension§!'@ 



E 



RK 



2m 
27 



J 2 Sx • S 2 



Si(x) - 



7T 



(44) 



where Si(x) is the integral-sine function 



Si(x) 



x sin(x) 



dx 



(45) 



and x = 2kpx (with x = \x2 — xi\) gives the phase of the characteristic oscillation of the 
interaction at twice the Fermi wavevector kp- 
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B. Green's functions in higher dimensions 



Once we have found the solutions for the purely ID Green's functions g v , which satisfy 
(p6|), we can find the solution of the planar extension of the equivalent problem, as demon- 
strated by ( P0|) , by means of Fourier transforms in the directions parallel to the barrier 

[(w - k\/2m) - H v (x)]g v (x, x'; kfu) = 5(x - x'). (46) 

Here again retarded and advanced versions of g^ 3D (x, x'; ky; uj) can be found simply through 
continuing analytically 

g r v /a (x,x';k\\;u) = g v (x, x'; k\\; u ± i8). (47) 

Correspondingly we find 

#<(x,x';k||;u;) = n F {u) [g*(x, x'\ kf uj) - g*(x, x'\ k y ; uj)] , rj G {L,R}. (48) 

It is important to realize at this point that for the retarded and advanced Green's functions 
g r l a the extension to higher dimensions only leads to a shift in the energy argument uj — > 
uj — H /2m, as can be seen from the form of their defining equation ([46]), i.e. 

g r r / a (x, x'; k f uj) = g^(x, x'; uj - Af /2m). (49) 

This property is seen to also translate in part to the functions g< with the only important 
difference that the energy arguments of the occupation functions remain unchanged. 

By using the properties (|48 ) and fl49|) one can establish that the Green's functions of the 



full system G r ^ and G^ can be represented in the following way: 

G [i)( x , x '\ k ll! = G r ^(x, x'\ uj - fcj/2m), (50) 
Gf Jx, x'\ ky; u;) = n^(o;)ri,(x, x'; uj — kj/2m) + rip(uj)T R (x, x'; uj — kj/2m), (51) 

where the functions T v , rj G {L, R} associated with the Fermi functions n F can be written in 
terms of sums and products of spatial derivatives of the g^ a , r] G {R, B, L} when (|30|)-(|37D 
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are used to divide G^(x,x') into contributions containing left or right occupation functions 
only. For example one finds: 

T L{R) (L, R) = -^1b(L, R) {j r %[R(L), R{L)] + Jb[R(L), R(L)}} 

{i a L {R) [L(R),L(R)]-YL {R) [L(R),HR)]}- (52) 



X 



The properties (^) and fl5"T[ ) prove to be useful for the evaluation of the frequency and 
wavevector integrations in the expression for the monolayer interaction (|2"2"|), which is also 
needed for the calculation of the slab interaction in (|23[) . If the interaction is considered at 
zero temperature and the additional property is used that the density of states factor coming 
from the 2D kn -integration in (|22"D is just a constant, one can reduce the three integrations 
in (|22|) to a sum of two single integrations as 



d 2 k» r°° duj 



(2tt) 2 J-oo 2tt 

2 m fli, fti-eV 



0(/i - w)f L (w - i) + e(A1 _ u _ e y)f - 



■ / dzifi — z)Tl(z) + / (izf/^ — 2 — eV)T R (z), (53) 
jo J- e y 

where the 0-functions derive from the sharp Fermi distributions n^uo) at T = and 

T v (z) = lm{V v (x 1 ,x 2 ; z) [G r (x 2 ,x 1 ; z) + G a (x 2 ,x 1 ; z)}} . (54) 

The lower limits at and — eV in ([53]) follow from the fact that there are no states below 
the bottoms of the conduction bands in the left and right lead, respectively. As a result 
of this simplification the interaction energy between monolayers of spins in 3D is not much 
harder to evaluate than the interaction energy in the purely ID case. 

For a 3D planar junction in equilibrium without a barrier we find from ([22] ) that the 
interaction density for monolayers is 



2mkp / 2D ^2 2D 
RK ~ 2(2^2 \ J Ps~d) ^ Si ' &2 



Si(x) - 



vr xcos(x) ~ sin(x) 
1 X 2 



(55) 



which is evidently similar to the purely ID result. The monolayer interaction is therefore 
often called quasi-one dimensionaS. The slab interaction from fl2"3] ) in this case assumes 
the formil 



E 



RK 



2m ( Jp 3 a E d 



v 



2D 



-Si • S 2 {F{s + 21) - 2F(s + I) + JF(s)} , 



(56) 



16(2vr) 2 

where = p^-dPs-d is the 3D density of s-ti spins within the slabs, s = 2d + R — L is 
the spacing between the slabs, I is the width of each slab (see Fig. |]) and T{x) is the range 
function 



Si(x) 



7T 



+ X cos(x) + sin(x) + 2Si(x), X = 2k F x. 



(57) 



Another special case worth noting in this context is the interaction between two magnetic 
half-spaces separated by a non-magnetic spacer. By letting I —>■ oo we find that (|5E|) simplifies 
to 



77ISOO 

-^RK 



2m(jp™ d y 



,2D 



16(2tt) 



Si-S 2 + 2 



Si(x s 



7T 



+ Xs cos(x s ) + sin(xs 



(58) 



where Xs = 2kps. 



In order to obtain solutions for (|T8| ), (p2|) and (^3[) also for more general cases we have 
to perform the corresponding energy integrals numerically. 



V. NUMERICAL RESULTS 

A. Comparative results in equilibrium 

Numerical results for the RKKY interaction across a tunneling junction in equilibrium 
have been given by Mukasa et a/.0 These authors were particularly interested in providing 
a theoretical model for possible applications in exchange force microscopy, as mentioned in 
the introduction, which is based on the RKKY interaction as the dominant force between 
a tunneling tip and a sample. For this purpose Mukasa et al. performed an approximate 
version of scattering wave perturbation theory for free electrons in a ID system by including 
the transmission coefficient of an electron tunneling through the barrier at an intermediate 
stage in their calculation. In contrast, our analytic expression for the RKKY interaction is 
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exact (to order J 2 ) since the constituent Green's functions include scattering by the barrier 
exactly. 



In order to compare our results in equilibrium with the ones in Ref. [T(| we have 
implemented our calculation with the same model parameters using a Fermi energy of 
fi = fi L = fi R = 5.0eV in equilibrium and a lattice constant a® = 2.50A, along with a Fermi 
wavevector kp< = 1.26 x 10 10 m _1 , implying a relative effective electron mass of m/m e = 1.20, 
where m e is the bare electron mass. Fig. |](a) shows plots representing the dimensionless 
interaction range function with 

*{x) = -7T {^J 2 Si " S 2 } 1 E RK (x), (59) 

where E RK (x) was obtained from (|T8|) as evaluated in (|38|). In Fig. 0(a), the impurities 
are considered to be fixed at the electrode-barrier interfaces, while the width of the barrier 
is continuously increased from 0.0A to 5. OA, to represent for example the height of an 
STM tip above a sample. Without the barrier (V /fi = 0.0), the range function reduces to 
$(x) = ty[ty/2 — Si(x)], where \ — 2kpX. 

We have also plotted in Fig. |3](a) the interaction for this case and for a barrier height 
of Vo/n = 0.5, where the system is in a scattering state. We find that as the barrier 
height is increased from zero, the interaction varies with a longer wavelength, corresponding 
to a decreasing relative wavevector between the top of the barrier and the Fermi level 



K F = y/2m(jM-V ). 

For Vo/fi > 1, the strength of the interaction decays exponentially with the width of 
the barrier, with an exponent that increases with Vo/fi, and the crossover into the antifer- 
romagnetic (AFM) regime is lost once Vo/fi > 1.3. In the interval of Vo/fi G [1.0, 1.3], the 
interaction still experiences one slight crossover to the AFM regime. This can be understood 
to arise from the nature of the transmission and reflection coefficients of the barrier, which 
are not purely exponential, but a mixture of hyperbolic functions. 

However, our results explicitly do not show the large oscillation of the interaction far 
in the ferromagnetic (FM) regime as was obtained in Ref. |10|, markedly for their curve 
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Vo/ fj, = 1.05. Such a behavior is unphysical in a genuine tunneling situation, and reflects the 
approximate treatment of the transmission of scattered waves through the barrier in Ref. [H]. 
Our Green's function approach in comparison includes the single particle barrier potential 
fully from the beginning and therefore allows an exact evaluation of the RKKY interaction. 

Other interesting cases to investigate are the interaction of two monolayers, or of two 
semi-infinite slabs of magnetic impurities in the presence of a tunneling barrier. Such systems 
could for example be realized by coating both sides of a fixed or mobile tunneling junction 
with a magnetic material. For these cases one can make use of (p2[) for monolayers and of 
(f2~2|) and (|23|) for slabs, both in conjunction with (|53|) . We shall consider (|2~2"D in equilibrium 



and set eV = in (|53|). In Fig. 0(b) and Fig. 0(c) we show results for the interaction range 
functions 

t (.Tn 2D ,Y u 2D ^ .Sj 

2(2tt) 



^) = --{^I(^) 2 ^Si-S4 (60) 



n 3D \ 2 ,,2D 



n 2m JpfJ v 2 
*?(*) = ~\ \ 6( J a S i- S 4 (61) 

for monolayers and for slabs, respectively, which were obtained for a planar version of the 
arrangement used for Fig. 0(a). From ([55]), for zero barrier height, $ m reduces to $ m = 
n{n/2 - Si(x) - [xcos(x) - Mx)}/x 2 } and $ s °° to $ s °° = -tt/2{[ X 2 + 2][Si( X ) - tt/2] + 
Xcos(x) + sin(x)}. The curves in Fig. 0(b) show that characteristically the monolayer 
interaction $ m decays much faster than the purely ID interaction $, shown in Fig. 0(a), 
when compared with equal parameter values. Especially for V /fi > 1.0, $ m shows a much 
stronger decay than $. This can be understood as a result of the fey-integration in ( |2"2] ) which 
effectively averages the interaction of one spin on one side of the junction with all other spins 
on the opposite side. This average which extends over many oscillatory contributions leads 
to destructive interference effects which results in the observed damping of $ m relative to 

When the RKKY interaction of infinite half-spaces represented by is compared to $ m 
and $ it is seen that the double spatial integration from (|23|) increases the relative oscillation 
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strength of when compared to $ m , and at Vq/[i = 1.05 the relative cross-over <3> s into 
the AFM regime becomes even stronger than the one for <3>. 

B. Nonequilibrium behavior 

1. Impurities on the electrode-barrier interfaces 

In the following we establish how the results shown in Fig. |^ for the interaction in 
equilibrium are modified when a finite bias is applied to the junction. Fig. |4](a) and Fig. 
^(b) show the ID interaction between two magnetic impurities placed on opposite electrode- 
barrier interfaces for various strengths of the bias in conjunction with initial equilibrium 
ratios Vo//i = 1.05 in Fig. ||(a) and Vo//i = 1.5 in Fig. |](b), and with otherwise the same 
model parameters as in Fig. |3|(a). 

As the bias is increased in Fig. |](a) the interaction starts to exhibit oscillations when the 
right edge of the barrier potential V(R) = Vq — eV is pulled below the chemical potential 
of the left hand side fi L . This oscillation arises because the wavefunctions of high energy 
electrons tunneling from the left exit the barrier through its sloping part, and become 
oscillatory over the distance where they are above the barrier. In this regime we find that 
the wavevector of the oscillation can be roughly approximated by the wavevector of electrons 
tunneling from the left Fermi level at the position of the right interface with the barrier, i.e. 

q F = ^2m[fi-V + eV}, (62) 

giving an oscillatory wavelength Ax = l-ftjqp. The wavelength of the interaction becomes 
smaller, i.e. the value of qp in (|62| ) increases, as the bias is increased. In Fig. f|(a), one can 
see that as the bias is turned up, the antiferromagnetic region initially vanishes and then 
almost reappears at a smaller distance between the impurities. 

In Fig. ^(b), for eV/fi = 1.50, a similar behavior to the one in Fig. |](a), for eV/fj, = 1.05, 
is observed. The main difference to Fig. ^(a) is that for eV/fi = 1.50 and in equilibrium 
the interaction just decays exponentially with no AFM region, whereas out of equilibrium 

22 



such an AFM region is established as the bias is turned up. The fact that the oscillations 
caused by high bias are not centered around the $ = line in Fig. |](a) and Fig. f|(b) can 
be understood to arise from the asymmetric shape of the sloping barrier. From these figures 
it is apparent that in an arrangement where the impurities are attached to the electrode- 
barrier interfaces, such as an STM tip and sample, both coated with magnetic materials, 
the interaction becomes tunable between FM and AFM by varying the bias alone. 

Since, as mentioned before, an interesting application of the present system would be 
in exchange force microscopy — where the force caused by the exchange interaction on a 
tunneling tip in an STM is measured — we have plotted in Fig. |4](c) and Fig. |](d) the spatial 
derivative —d$>(x)/dx of the range function from Fig. f|(a) and Fig. f|(b), respectively. 
Both Fig. |](c) and Fig. |](d) show explicitly that the onset of oscillations in the interaction 
for an appropriate bias will lead to a force which is alternating in sign. An estimate of the 
absolute strength of the exchange force and predictions that it should be measurable with 
a state of the art STM are postponed to a discussion in Sec. VI. 

We next consider the effect of a finite bias on a system of two magnetic monolayers 
interacting across a 3D planar barrier. In Fig. |](e) and Fig. |](f) we plot the range function 
$ m , using fl2"2]), Q53]), and (pD|), for a planar version of the arrangement used in Fig. [|(a) 
and Fig. §(b). As in Fig. |](a) and Fig. |](b) the interaction starts to exhibit an oscillatory 
behavior once the slope of the barrier gets steep enough. In both Fig. f|(e) and Fig. f|(f) it 
is evident that the oscillations can reach into the AFM region, although, as in Fig. f|(a) and 
Fig. ||(b), the oscillations are again not centered about the $ m = line, but are shifted into 
the FM region. 

One more remark is in order when nonequilibrium results for monolayers are compared to 
results in equilibrium. In equilibrium the total interaction energy E^ K in ( f22"|) is proportional 
to /i, and we have normalized the range function $ m in Fig. [|(e) and Fig. |](f) with respect 



to the equilibrium Fermi energy \i. However, the total nonequilibrium interaction from ( |22| ) 
depends in magnitude on a non-separable mixture of /z, eV, /i + eV and /i — eV. This would 
be evident if we had plotted the interaction back to x = where the results for various 
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strengths of eV would no longer converge in a single point, as they do in equilibrium for 
different Vo//x. 

Altogether, our results in Fig. |](a) and Fig. ^(b) as well as in Fig. |](e) and Fig. [|(f) 
show that when the impurities are attached to the electrode-barrier interfaces a switching 
of the interaction can be achieved in many situations by changing the bias alone, but in 
general this switching behavior depends quite strongly on the particular properties of the 
barrier. In the following we will show that such a switching is much more reliably achieved 
by placing the single spins or layers of spins within the electrodes a finite distance away 
from the interfaces with the barrier. 

2. Impurities within the electrodes 

We now consider the nonequilibrium behavior of the interaction when the impurities are 
placed inside the electrodes a distance d away from the electrode-barrier interfaces. In Fig. 
|5](a) we show a surface plot of a ID arrangement with the relative barrier height Vq/ fi = 1.5, 
where the barrier width is kept constant at 1.0A [cf. corresponding point in Fig. |](b)]. The 
impurities are now moved away from their initial positions on the electrode-barrier interfaces 
to a maximum distance of d = 5. OA from either interface (plotted across the figure). At the 
same time the bias is increased from eV/fi = 0.0 to eV/p, = 1.5 (plotted into the depth). We 
have overlayed a contour plot to make it easier to identify which regions of the surface lie in 
the FM or AFM regime. The solid zero contour line indicates the boundaries of these regions. 
In the same way, we show in Fig. |5|(b) the interaction density between two monolayers and 
in Fig. H(c) the interaction density between two slabs of finite thickness which are placed a 
finite distance inside two 3D electrodes. The range function in Fig. |5|(b) is again normalized 
with respect to the equilibrium Fermi energy fi. The case of interacting slabs turns out to 
have quite special features which will be discussed later in this section. 

One can see that in both Fig. |5|(a) and Fig. |^(b) the interaction is already oscillatory in 
equilibrium as d is varied, even though V / fi > 1.0, since now the electrons have to travel over 
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a finite region on either side of the junction where their wavefunctions are oscillatory. The 
presence of the barrier in this case leads to an overall exponential damping which reduces 
the strength of the oscillations everywhere. As the bias is turned up these oscillations evolve 
into an interference between up to five contributing components which can be explained 
as follows: To the order of perturbation theory considered, the spins on the left side of 
the junction interact with the ones on the right through electrons tunneling between the 
locations of the spins. Electrons in the vicinity of the left spins which are able to perform 
this process are available up to the Fermi level on the left. The wavevector ki,L of the 
spin polarization of the conduction electrons in the left lead is determined by the cut-offs 
of the frequency integration in ( |i8|) and (|22"D at the Fermi energy and at the band bottom 



on the left side V , so that kiL = J2m(/j, L — V L ). This wavevector is, however, a different 



one, namely km = y 2m(/z L — V R ) once the tunneling electrons have penetrated the barrier 
and interact with the spin on the right side. The same process applies to the spins on the 
left feeling the presence of the ones on the right. The corresponding wavevectors involved 



therefore comprise kij = y2m(/i* — Vi), G {L,R}), i.e. four different ones in principle. 
In addition to this, one also observes a quite strong oscillation which contains the wavevector 



k e y = v2rneV. This oscillation can be understood to arise from the interval in the frequency 
integration in the expression for the interaction (18 ) which extends from the band bottom 



of the right lead to the band bottom of the left lead, i.e. over V L — V R = eV. Since all cases 
of the interaction studied here have a one- dimensional or quasi-one-dimensional behavior, 
both the densities of states of the electrons in the left lead and in the right lead exhibit a 
singular behavior at the respective band bottoms. The existence of the singular parts of 



these densities of states at the integration limits gives rise to the y/2meV -oscillation, which 
in some situations becomes the principal oscillation in the nonequilibrium system. 

Assuming that /jl r moves by the same amount eV as V R , and with V L normalized 
to zero as shown in Fig. [I], this set of wavevectors reduces to a total of four, namely 
fee/ y / 2m / u L , y / 2m(/i L + eV), ^2m(fi L — eV)j due to the sharp Fermi distributions n L 
and n R , and the k e y = \/2meV oscillation as discussed before. In both Fig. |5](a) and Fig. 
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H(b), a superposition of these principal wavevectors can be found. Once the bias eV moves 
the Fermi energy on the right below the band bottom on the left, fi R < 0, the corresponding 
components of the oscillations turn into an exponential decay, leaving only three contribu- 
tions to the oscillations. This transition can be seen as a kink in the contour plots in Fig. 
|5](a) and Fig. |5](b) at a bias of eV = 1.0. At higher biases, for arrangements where the 
impurities lie very close to the barrier, i.e. when d/(R — L) <C 1.0, barrier effects such as 
those displayed in the plots of Fig. |](b) an d Fig. ||(f) become more important. However, for 
ratios d/(R — L) > 1.0 we expect these effects to be minimal when the width of the barrier 
is held constant. Both Fig. |5|(a) and Fig. |5|(b) show explicitly how the interaction between 
two single spins in ID and between two layers of spins, respectively, is tunable between 
ferromagnetic and antiferromagnetic coupling at a fixed impurity configuration by varying 
the bias alone. Furthermore it is clear in Fig. |5](a) and especially in Fig. ^(b) that the in- 
teraction falls off significantly less rapidly in the nonequilibrium regime. This circumstance 
has particularly strong consequences for the interaction between slabs of spins. 

In Fig. |5|(c) we have plotted the interaction between two finite magnetic slabs of spins 
with a width of the slabs of I = 10A using otherwise the same conditions as were taken for 
the monolayer-case shown in Fig. |5|(b). While it can be seen in Fig. |3](c) that in equilibrium 
the interaction for slabs converges to a finite value in the limit when the thickness of the slabs 
goes to infinity, I —>■ oo, this is no longer the case out of equilibrium. Rather the interaction 
is seen to become roughly proportional to the thickness of the slabs for thicknesses I ^> ir/hp- 
This can be attributed to the longer range of the monolayer interaction as mentioned before, 
which destroys the quite sensitive convergence of the double integral in (^) for I — > oo. 

When the slab interaction (|23|) is calculated one can exchange the spatial integrations 
with the frequency and wavevector integrations. This has the advantage that the spatial 
integration can be performed analytically first. Since the explicit representation of the 
integrand is quite involved, we have postponed it to Appendix A. From ( |A10| ) it can be seen 
that the integrand contains several terms which are proportional to /. In equilibrium, as 
shown in (TA"8|), these terms cancel one another, but out of equilibrium, when the frequency 
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integration for these terms is cut off at different points through different Fermi occupation 
functions, this cancellation ceases to be complete. Therefore, when the bias is switched on, a 
residual /-dependence remains, which increases as the bias is increased. The total interaction 
energy for slabs therefore scales extensively with the width of the slabs. It should be noted 
that this extensive dependence of the slab interaction energy out of equilibrium is not caused 
by the offset of the band bottoms in our model. If an equilibrium system is considered in 
which materials with different band structures form a junction, the slab interaction continues 
to converge always [cf. (|A8|)]. 

For this reason we have compressed the scale of the plot in Fig. |5|(c) in the direction of 
the bias by a factor of [1 + 40 x eV/fi}' 1 , so that, as shown, the height of the oscillations 
stays approximately the same. Moreover, one can see that the kpd oscillations observed 
in equilibrium vanish completely at a very small bias and turn into oscillations with phase 



(d/2)V2meV. That this is the case can be seen from the overlayed contour lines, which, 



particularly in the upper right corner, show a dependence oc 1/ y / 2meV when taken as a 
function of eV. Mathematically, the dominance of these oscillations can be understood 
to arise from the spatial integrations in formula (^). When the spatial integrations are 
interchanged with the frequency integration in (^), as it is done in ( |A1| ) in Appendix A, the 
space integrations are effectively taken over products of the oscillatory, i.e. trigonometric, 
wavefunct ion- like expressions Vl{r) from (|A^) which relate to electrons in the left and right 
leads [cf. ( |A4|) -(|A7j)l. This yields an extra factor of (i?^))" 1 i n Cl(r) from ( |A5|) . A factor 



of Cl(r) occurs at least once in every term in the integrand of the subsequent frequency 
integration in (|A10|) , which enhances the peak in the quasi-lD density of states of the 



monolayer system occurring on the bottom of the band in the corresponding lead. The 
frequency integration itself very much acts as a Fourier transform into real space, which gives 
the peaks the effect of frequency components of the real-space oscillations of the interaction. 
These peaks are just an energy interval eV apart, which explains the occurrence of the 
k e v = \f2meV wavevector as a difference wavevector between these components. The 
strength of k e y is then determined by the strength of the peaks, i.e. the amplitude of the 
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Fourier components. 

Since, as is shown in Fig. [|(c), the interaction between magnetic slabs oscillates with 
almost only the k e y contribution present, it should be controllable in a simple fashion by 
tuning the bias. 



VI. DISCUSSION 

In this section, the possibilities for the experimental study of the RKKY interaction in 
tunneling systems out of equilibrium are discussed. As mentioned in the introduction, one 
possible way to probe the behavior of the interaction that we predict is to use an STM to 
measure the exchange force, such as shown in Fig. |](c) and Fig. §(d). 

In order to estimate the absolute value of the exchange force, we must estimate J in the 
prefactor of the force function between two magnetic impurities, adsorbed to the STM tip 



and the sample, respectively. One method to estimate J, which is also used in Ref. [L0[ is to 
use the equation for the Kondo temperature T% 

1 



T K = — exp 
k B 



JgZM 



(63) 



where Kb is the Boltzmann constant and gl® d ([i) = (27r) _1 y 2m//i is the ID density of s-d 
spin states at the equilibrium Fermi energy fi. For measured TVs of about Tk ~ 100 K 
one obtains Jgl® d (fi) ~ 0.2. From this the function F(x) representing the total measured 
exchange force is determined as 

From the parameters introduced above, we find 

F ~ Si ■ S 2 ^^ x O (lO- 20 Nm) , (65) 

which for the range functions plotted in Fig. |](c) and Fig. ^(d) leads to forces of about F ~ 
(9(10 -11 N). As noted in the literature0il the resolution of current atomic force microscopy 
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is about 10 _11 N — 10~ 13 N, which means that the exchange force should be experimentally 
observable. 

For the observation of the RKKY interaction in layered tunneling systems which exhibit 
giant magnetoresistance phenomena the results shown in Fig. [| should give an appropriate 
prediction. Especially for the case of 3D systems the results for the interaction between 
interacting slabs of impurities shown in Fig. H(c) would be most applicable. The most 
surprising result in this situation is that we find that the interaction scales extensively with 
the thickness of the slabs. In principle the expression for the interaction we study would 
diverge in the limit of semi-infinite slabs, but of course this would not be observable in real 
systems, since processes like spin relaxation or the scattering of electrons on non-magnetic 
impurities would eventually impose a maximum range on the interaction. However, for slabs 
of reasonable thickness the switching of the bias should produce an observable change in the 
decay behavior of the interaction, i.e. in a system with fixed spin positions an anomalous 
increase of the interaction strength should be measurable. 

Furthermore, the tunability of the RKKY interaction in a tunneling system out of equi- 
librium is particularly interesting in view of the fact that the impurity configuration is 
normally pre-set in a solid system, which means it is usually not possible to directly observe 
the oscillatory dependence of the interaction on the impurity spacing. With the present 
arrangements it becomes possible to observe oscillations of the RKKY coupling via the vari- 
ation of wavevectors on either side of the junction as the bias is varied. If one follows for 
example the direction of the bias at a given distance of the spins in any of the contour plots 
shown in Fig. |5], one can almost always observe at least one crossing from FM to AFM 
coupling or vice versa. 

Since, out of equilibrium, as indicated before, the interaction depends on several con- 
tributing wavevectors, the relative dependence on these could for example be investigated 
experimentally by using different materials for the connecting leads and by placing the spins 
asymmetrically around the barrier, e.g. on the electrode-barrier interface on one side of the 
junction and further within the electrode on the other side. As the arrangements in Fig. 
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H(b) and Fig. ||(c) describe routinely achievable physical systems, this suggests that the 
phenomena found theoretically for such systems should also be accessible to experiment. 

VII. CONCLUSION 

This paper presents for the first time a theoretical treatment of the RKKY interaction 
in systems out of equilibrium. Except for processes involving ultrashort time dependent 
excitations, a proper nonequilibrium situation seems only routinely attainable in structured 
systems which include a potential barrier. The main achievement of this work is that we 
have obtained a proper field theoretic description of such a system which is much more 
systematic than conventional scattering wave perturbation approaches to the problem (cf. 
Ref. Difficulties such as taking proper account of the different occupation functions 
in different parts of the system are overcome as well as the problem of how to normalize 
the wavefunctions involved, since the Green's functions used in the present description are 
always properly normalized. Our treatment is adaptable to the inclusion of further many- 
body effects in the problem, such as carrier-carrier interactions in the electrodes and the 
interaction with phonon modes. 

One effect of biasing the system out of equilibrium is that the oscillatory exchange in- 
teraction in various dimensions exhibits strong interference effects, leading to one or more 
changes of the type of coupling (between FM and AFM) at a given impurity configuration 
as the bias is varied. This behavior arises as a result of an interference between several 
fundamental oscillations due to a mixing of different wavevectors. The possibility of tuning 
the interaction through changing the bias alone could become an important effect in ap- 
plications of nonlinear switching devices using layered magnetic structures with a potential 
barrier. Another important effect is that out of equilibrium the range of the interaction 
increases, leading to an interaction energy that scales extensively with the system size in the 
direction perpendicular to the barrier for interacting slabs of spins. A closer study of this 
phenomenon for the 3D case presented here, as well as for the ID case of interacting lines 
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of spins, where the interaction in equilibrium is known to lead to a helical ordering within 
each line, is being undertaken. 

Our results should be applicable to a very broad variety of conceivable structures and the 
formalism we have presented here is particularly suitable to be adapted to such situations. 
Extensions to include varying effective masses and other material properties such as band 
structure and different band filling in the various sub-parts of the system are obvious. A 
particular example for a possible extension would be to study the interaction in double 
or multi-barrier systems out of equilibrium. Furthermore, in order to facilitate a direct 
comparison to experimental results a next step could include the calculation of the RKKY- 
perturbed spin-polarized tunneling current across systems of this kind. We also hope our 
work will encourage the experimental study of giant magnetoresistance phenomena out of 
equilibrium where one can equally expect interesting interference phenomena to occur as a 
result of the different relative distances of the Fermi surfaces to the bottoms of the conduction 
bands involved. 
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APPENDIX A: NONEQUILIBRIUM BEHAVIOR OF THE INTERACTION 

BETWEEN SLABS 



In the expression for the RKKY interaction between 3D slabs of spins from (|23|) , we 
exchange the spatial integrations with the frequency integral, to obtain 

(jpf d )V D ^ doJ dk 2 L _ d fR+(d+l) 

e rk = V ,„ \ A S x ■ S 2 / — / / _ dxj dx 2 (Al) 



(2m) 4 J-oo 2tt J (2tt) 2 Jh-(d+i) Jn+d 
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Im 



< r 



r lL (x 1 )G (Si) (L 1 R)r, R (x2)\ [vU^G^R, L)rf L ( Xl ) + Vr (x 2 )G{ 0) (R, L) V l{x x ) 

In ( |A1[ ) we have used (|5^ ) and Q55D to rewrite the expression G^(xi, £2) 
G ? (q- ) (x2, xi) + G" - ) (x2,xi) in the integrand. In addition, we have introduced 



r/a 



X0=L{R) = d ^tR)^X l{2) ] 
= + (_)2mexp {T(±W L { R) [x m - L(R)]} , (A2) 

where g r ^ R ^(x,x') is taken from (|39|). Furthermore, the explicit frequency and wavevec- 
tor dependence of the integrand in ( |A1| ) was omitted for brevity. The expression for 
G^(xi, X2) G( )(x2, x\) + G^(x2, xi) in the integrand can now be further grouped as fol- 
lows: 



x =L(R) 



(A3) 



Gf 0) (x 1 ,x 2 ) [G( 0) (x 2 ,xi) + G" 0) (x 2 ,a;i) = 
(2m)- 4 {G5 ) (L, J R)G[ 0) ( J R,L)[r ? l(x 1 )] 2 ^(x 2 )^(x2) 

+ G< ) (L, J R)q 0) ( J R,L)^(x 1 )r ? 2(x 1 )[^(x 2 )] 2 

+ G[ 0) (L, i2)q 0) (i?, L) [rf L {x l )fr ]R {x 2 )r 1 < (x 2 ) 

+ q 0) (L, J R)q 0) (i?,L)r ? <(x 1 )r / 2(a:i)[^(x2)] 2 

+ q 0) (L, J R)G^ 0) ( J R,L)r / 2(x 1 )r / 2(x 1 )^(x2)r ? <(x 2 ) 

+ Gf 0) (L,i?)q 0) (i?, J L)r 7 <(x 1 )^(x 1 )^(x 2 )^(x 2 )}. 

When now the functions are replaced by their definitions in terms of retarded and ad- 



vanced functions, 



Vl(r) 



n 



L(R) 



Vl(r) ~~ ^Il(r) 

integrations can be accounted for by introducing the following terms: 

rL-d / pR+d+l 



one can see that all occurring dxi and dx 2 



,L(R) 



L-d-l \JR+d 

(2m) 2 i 



dx l(2)V r L(R)( X l(2))V a L(R)( X l(2)) 



(A4) 



1l(r) Ql{r) 

(2m) 2 



{exp [i(q r L{R) ~ q a L{R) )(d + l) 



— exp 
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Il(r) 



L-d 



R+d+l 



L-d-l \JR+d 

(2m) 2 i r 
"F"~w^ i ex P 



2n r,a 



dx l(2)[V r L {R)(Xl(2))] 2 
±2i QL(R)( d + l )] ~ eX P 



(A5) 



±2iq r L / ( a R) d]} 



While from ( |X4| ) produces a term proportional to I for energies above the band bottom 
in the respective lead, it is seen from ( |A4| ) and (JA3|) that 



A) 



_ tr/a 
L(R) — ?L(R)' 



V L(R) < o. 



(A6) 



From here we find 



L-d rR+(d+l) f 

dxi / dx 2 Gf 0) (xi,a;2) G r (Q) (x 2 , x x ) + G a m {x 2 , x x ) 

L—(d+l) J R+d 



(A7) 



(2m)- 4 {Gf 0) (L, i?)q )(^ ^)S€r + G fo)(4 #)<^o)(^ 
+ q 0) (L, ^)q 0) (^ [& - &] + G a {Q) (L, R)G a {Q) (R, L)n L F (u)& - & 

+ G{ 0) (L, R)G a (0) (R, L)nf H£ " &] + G a {0) (L, R)G[ 0) (R, L)n L F {u)^ [£ - &]} . 

Since the problem we consider has time- reversal symmetry, G r ^(L,R) = G r ^(R,L) al- 



ways holds. When (|A7|) is considered in equilibrium, we know that also G^ (L, R) 



n F (u ) 
to 



G to) (A R) - Gr (o) i L i R) 



(0)> 



holds. In such a situation it is seen that 



simplifies 



(A8) 



L-d i-R+(d+l) f 

dxi dx 2 Gf 0) (x 1 ,x 2 ) G\ Q) (x 2 ,xi) + G a {0) (x 2 ,X] 

n F (u) { [Gf 0) (L, R)} 2 ClCr - [q 0) (A R)] 2 



where it was assumed that £l and £r may also be different in equilibrium, e.g. if the materials 
on the left and right of the junction have different work functions. Clearly in (|A8|) all factors 
proportional to I have vanished as expected. When I is taken to / — > oo the functions Cl(r) 
oscillate with I, as is seen from ( |A5| ) , and eventually reduce to 



(2m) 2 i 
2n r/a 

Z( iL(R) 



exp 



±2iq r L / ( a R) d 



(A9) 



due to the retarded/advanced property of the q^m^, which causes an exponential decay in 
the relevant functions at large distances. 
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Out of equilibrium, however, the cancellation of the /-proportionality ceases to be com- 
plete and the expression for the slab interaction (|23| ) can be rewritten by means of (5^) and 
(|53|) for zero temperature as 



E RK 



(Jp[ 



l7r 2 (2my 



Si-S 5 



(A10) 



x <! J dz(fi — z)lm 

dz(fi — z — e\^)Im 

eV 



where we have introduced 



(T L (z) - Gfa(zj) E(z) + (Gfo)W) 2 ''&(z)C R (z 

(T R (z) + G[ Q) (z)) E(z) - (G r {0) (z)f &(z)C R (z) 



E(z) = G\ Q) {z)H{z)e R {z) + Gfo(z)&(z)#(z) 



(All) 



with the abbreviation: G(q\ (L, R; z) = G(p)(z) and r L ( R \(L, R; z) = T L ^ R )(z). Note that 



in ( |A10| ) terms containing C,l^rGIq){z)G^(z) have vanished, since they are entirely real. 



From the definition of the ^I(R) m ( 1^1 ) ^ * s seen that the term H(z) from ( |A11| ) is pro- 
portional to I for z > 0. When — eV < z < 0, H(z) always has one /-proportional and 
one exponentially decaying component, where the latter describes the penetration of low- 
lying electron states in the right lead to the position of the spins in the left one. The 
/-proportional component, however, vanishes for energies below zero, since in this case the 
term T R (z) + G r , Q \(z) G r ^(z)^ 7 L (z)C, R (z) turns out to be entirely real. The remaining /- 
dependence of subsequently translates to the slab interaction from ( A10|) for energies 
z > 0, where it is seen to produce terms proportional to /, i.e. which scale extensively with 
the thickness of the slabs, as well as the equilibrium terms from (KB) before, which show no 
/-proportionality. 
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FIGURES 

FIG. 1. Schematic representation of the one dimensional tunneling system described in the 
text. Two magnetic s-d impurities are situated in the left and right lead at equal distances d away 
from the electrode-barrier interfaces L and R, respectively, and interact with each other through 
the tunneling of electrons across the barrier. The points L and R at the same time serve as the 
partitioning points within the Caroli/Feuchtwang formalism, indicated by the two dashed vertical 
lines. 

FIG. 2. Cross-sectional drawing of a magnetic multilayer structure, obtained as a 3D planar 
extension of the ID system shown in Fig. 1. Two magnetic slabs M are separated by non- magnetic 
spacer layers S and a planar tunneling barrier B extending between points L and R. The spins in 
the magnetic slabs are assumed to have the same orientation within each slab due to a dominance 
of ferromagnetic (FM) coupling at short distances. 

FIG. 3. Range functions: (a) for interacting magnetic impurities in ID from (59), (b) 

& m (x) for interacting magnetic monolayers in 3D from (60), and (c) <E^°(x) for semi-infinite in- 
teracting magnetic slabs in 3D from (61) plotted against the distance x = R — L between the 
spins in equilibrium, where the spins are considered to be fixed on either interface of the barrier 
(d = 0). The barrier width is increased from O.oA to 5.0A and the height of the barrier is increased 
throughout the plots as indicated. For Vo/fj, = 0.0 (a), (b) and (c) show the range functions <&(x), 
$ m (x) ^^°(x), respectively, for a free electron system. 

FIG. 4. (a),(b): range function $(x) from (59), (c),(d): force function —d&(x)/dx, and (e),(f): 
range function $ m (x) from (60) for finite bias eV . The initial relative barrier height in equilibrium 
is fixed at (a),(c),(e): Vo/fx = 1.05 and (b),(d),(f): Vo/fx = 1.50 with otherwise the same system 
parameters as in Fig. 3 (fx is the Fermi energy of the leads in equilibrium). In four steps a bias eV 
of up to eV/n = 2.0 is applied to the junction. 
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FIG. 5. Surface plot of the range function (a): (b): $ m and (c): & s = [1+40 x eV / fi}- 1 &- 10 
when the impurities are moved within the electrodes. In the present arrangement the barrier height 
is fixed to Vo/n = 1.50 and the barrier width to R — L = 1.0 A. The distance d of the impurities 
on either side of the barrier is increased from d = (0.0 — 5.0) A (plotted across), while at the same 
time the bias is varied from eV = (0.0 — 1.5)eV (plotted into depth). The thickness of the slabs in 
(c) was taken to be I = 10 A. In (a), (b) and (c) the front faces of the surfaces show the equilibrium 
range functions, which are oscillatory as the spins are moved through the leads. In order to be 
able to better identify the boundaries between FM and AFM regions, we have overlayed a contour 
plot showing the zero coupling (<E> = 0) contour. 
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Figure 4 
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